/**
  ETFAna project, Anyang Normal University && IMP-CAS
  \class ETFUserAnaTask
  \brief user class for tasks in the tracking stage
  \author SUN Yazhou, asia.rabbit@163.com
  \since 2024-01-06
  \date 2024-01-06 last modified
  \attention
  changelog
  <table>
  <tr>  <th>Date         <th>Author      <th>Description                   </tr>
  <tr>  <td>2024-01-06   <td>Asia Sun    <td>file created                  </tr>
  </table>

  \copyright Copyright (c) 2021-2024 Anyang Normal U. && IMP-CAS with LGPLv3 LICENSE
*/

#include <algorithm>
#include <TFile.h>
#include <TTree.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TClonesArray.h>
#include "ETFPXISection.h"
#include "ETFUserAnaTask.h"
#include "ETFDetectorManager.h"
#include "ETFTOFWall.h"
#include "ETFDCTa.h"
#include "ETFDC.h"
#include "ETFPlaHit.h"
#include "ETFDCHit.h"
#include "ETFHPTDCCh.h"
#include "ETFAmpCh.h"
#include "ETFMsg.h"
#include "ETFMUSIC.h"
#include "ETFTimer2.h"
#include "ETFTimer4.h"
#include "ETFT0Array.h"
#include "ETFVMEEvent.h"
#include "ETFPXISection.h"
#include "ETFTrackManager.h"
#include "ETFTaTracker.h"
#include "ETFPostMagTracker.h"
#include "ETFTrack.h"
#include "ETFAnalyzer.h"
#include "ETFTrackTa.h"
#include "ETF3DTrack.h"
#include "ETFPID.h"
#include "ETFMath.h"

using std::swap;
#define er ETFMsg::Error
#define within ETFMath::Within
#define lii ETFMath::lii

ETFUserAnaTask::ETFUserAnaTask(const char *name) : ETFTask(name), fT1p(0), fT1v(0),
fT0p(0), fT0v(0), fT0F1(0), fTOFW(0), fT0Arr(0), fMU{0}, fdcTa{{0}}, fdc{0}{
  et = new TTree("et", "Extra Terrestrial");
  fBrho0 = 4.9896745; // Ni58 setting, 30Dec, 2023
} // end ctor

ETFUserAnaTask::~ETFUserAnaTask(){}

//////////////////////   Let's create some histograms   //////////////////////
void ETFUserAnaTask::CreateHistos(){
  Configure(); // get some detectors
  // configure the tree //
  // delta-E related
  et->Branch("dE", dE, "dE[2]/D");
  et->Branch("z0", &z0, "z0/D");
  et->Branch("z1", &z1, "z1/D");
  et->Branch("pileUpN", &fPileUpN, "pileUpN/I");
  et->Branch("pileUpT", &fPileUpT, "pileUpT/D");
  // tracking related
  et->Branch("k0", k0, "k0[2]/D");
  et->Branch("b0", b0, "b0[2]/D");
  et->Branch("k1", k1, "k1[2]/D");
  et->Branch("b1", b1, "b1[2]/D");
  et->Branch("k2", k2, "k2[2]/D");
  et->Branch("b2", b2, "b2[2]/D");
  et->Branch("pT1", pT1, "pT1[2]/D");
  et->Branch("pTAU", pTAU, "pTAU[2]/D");
  et->Branch("pTAD", pTAD, "pTAD[2]/D");
  et->Branch("pdcTaU", pdcTaU, "pdcTaU[2][2]/D");
  et->Branch("pdcTaD", pdcTaD, "pdcTaD[2][2]/D");
  et->Branch("pdcTaUWire", pdcTaWire[0], "pdcTaUWire[2][2]/D");
  et->Branch("pdcTaDWire", pdcTaWire[1], "pdcTaDWire[2][2]/D");
  et->Branch("pdc", pdc, "pdc[3][2]/D");
  et->Branch("ptofw", ptofw, "ptofw[2]/D");
  // f1 related
  et->Branch("x_f1", &x_f1, "x_f1/D");
  et->Branch("t0_f1", &t0_f1, "t0_f1/D");
  et->Branch("t0arr_m", &t0arr_m, "t0arr_m/D");
  // post-TA pid related
  et->Branch("tof1", &tof1, "tof1/D");
  et->Branch("tof1_v", &tof1_v, "tof1_v/D");
  et->Branch("aoz0", &aoz0, "aoz0/D"); // aoz in ribll2
  et->Branch("tof2", &tof2, "tof2/D");
  et->Branch("rho", &rho, "rho/D");
  et->Branch("beta0", &beta0, "beta0/D");
  et->Branch("bg0", &bg0, "bg0/D");
  et->Branch("brho0", &brho0, "brho0/D");
  et->Branch("beta", &beta, "beta/D");
  et->Branch("bg", &bg, "bg/D");
  et->Branch("aoz", &aoz, "aoz/D");
  et->Branch("aoz0c", &aoz0c, "aoz0c/D");
  et->Branch("z0c", &z0c, "z0c/D");

  // create some histos //
  if(!(fFolder = dynamic_cast<TDirectory *>(fRootFile->FindObject("UserAnaTask"))))
    fFolder = fRootFile->mkdir("UserAnaTask", "Statistics for the Analysis Stage");
  if(!fFolder) ETFMsg::Error("ETFUserAnaTask", "CreateHistos: fFolder nullptr");

  hT0DtPvV = new TH2F("hT0DtPvV", "Tstart L-R Time to Trig (VME vs PXI)",
    300, -5., 5., 300, -5., 5.);
  hT0DtPvV->GetXaxis()->SetTitle("PXI T0 tl-tr [ns]");
  hT0DtPvV->GetYaxis()->SetTitle("VME T0 tl-tr [ns]");
  hT1DtPvV = new TH2F("hT1DtPvV", "Tstop L-R Time to Trig (VME vs PXI)",
    300, -5., 5., 300, -5., 5.);
  hT1DtPvV->GetXaxis()->SetTitle("PXI T1 tl-tr [ns]");
  hT1DtPvV->GetYaxis()->SetTitle("VME T1 tl-tr [ns]");
  ht0arrMulti = new TH1F("ht0arrMulti", "t0 array multiplicity with fired status == 4",
    53, -1.5, 51.5);
  htofwMulti = new TH1F("htofwMulti", "TOF Wall multiplicity with fired status == 4",
    18, -1.5, 16.5);
  // to check the status of mtdc pile-up solution //
  hT0LR = new TH2F("hT0LR", "Tstart L : R - for checking mtdc pileup",
    300, 0, 47000, 300, 0, 47000);
  hT1LR = new TH2F("hT1LR", "Tstop L : R - for checking mtdc pileup",
    300, 0, 47000, 300, 0, 47000);
  hT0LE = new TH2F("hT0LE", "Tstart L : # of edges - for checking mtdc pileup",
    10, -1.5, 8.5, 300, 0, 47000);
  hT0RE = new TH2F("hT0RE", "Tstart R : # of edges - for checking mtdc pileup",
    10, -1.5, 8.5, 300, 0, 47000);
  hT1LE = new TH2F("hT1LE", "Tstop L : # of edges - for checking mtdc pileup",
    10, -1.5, 8.5, 300, 0, 47000);
  hT1RE = new TH2F("hT1RE", "Tstop R : # of edges - for checking mtdc pileup",
    10, -1.5, 8.5, 300, 0, 47000);
  hT0L = new TH1F("hT0L", "Tstart L of all Edges - for checking mtdc pileup",
    1000, 0, 47000);
  hT0R = new TH1F("hT0R", "Tstart R of all Edges - for checking mtdc pileup",
    1000, 0, 47000);
  hT1L = new TH1F("hT1L", "Tstop L of all Edges - for checking mtdc pileup",
    1000, 0, 47000);
  hT1R = new TH1F("hT1R", "Tstop R of all Edges - for checking mtdc pileup",
    1000, 0, 47000);
  htof1 = new TH1F("htof1", "TOF: T0(tstart) -> T1 (tstop)", 1000, 110., 180.);
  htof1->GetXaxis()->SetTitle("tof1 [ns]");
  htof1_v = new TH1F("htof1_v", "TOF: T0(tstart) -> T1 (tstop) from VME", 1000, 110., 180.);
  htof1->GetXaxis()->SetTitle("tof1 (VME) [ns]");
  htof1p_v = new TH2F("htof1p_v", "TOF: T0(tstart) -> T1 (tstop) PXI vs VME", 300, 110., 180., 300, 110., 180.);
  htof1p_v->GetXaxis()->SetTitle("tof (PXI) [ns]");
  htof1p_v->GetYaxis()->SetTitle("tof (VME) [ns]");
  htof2 = new TH2F("htof2", "TOF: T1(tstop) -> TOF Wall", 19, -1.5, 17.5, 500, 30., 120.);
  htof2->GetXaxis()->SetTitle("strip Id"); htof2->GetYaxis()->SetTitle("tof2 [ns]");
  hdE0 = new TH1F("hdE0", "music 0", 1000, 0., 4.3);
  hdE1 = new TH1F("hdE1", "music 1", 1000, 0., 4.3);

  htof1_v_dE0 = new TH2F("htof1_v_dE0", "tof1(VME) - MUSIC0", 300, 124., 153., 300, 0., 4.);
  htof1_v_dE0->GetXaxis()->SetTitle("tof1 from VME [ns]"); htof1_v_dE0->GetYaxis()->SetTitle("MUSIC0");
  htof1_v_dE1 = new TH2F("htof1_v_dE1", "tof1(VME) - MUSIC1", 300, 124., 153., 300, 0., 4.);
  htof1_v_dE1->GetXaxis()->SetTitle("tof1 from VME [ns]"); htof1_v_dE1->GetYaxis()->SetTitle("MUSIC1");
  htof1_dE0 = new TH2F("htof1_dE0", "tof1(PXI) - MUSIC0", 300, 124., 153., 300, 0., 4.);
  htof1_dE0->GetXaxis()->SetTitle("tof1 from PXI [ns]"); htof1_dE0->GetYaxis()->SetTitle("MUSIC0");
  htof1_dE1 = new TH2F("htof1_dE1", "tof1(PXI) - MUSIC1", 300, 124., 153., 300, 0., 4.);
  htof1_dE1->GetXaxis()->SetTitle("tof1 from PXI [ns]"); htof1_dE1->GetYaxis()->SetTitle("MUSIC1");

  htof2_dE1 = new TH2F("htof2_dE1", "tof2 - de1", 500, 30., 80., 500, 0., 4.3);
  htof2_dE1->GetXaxis()->SetTitle("tof2 [ns]"); htof2_dE1->GetYaxis()->SetTitle("dE");

  // t0arr statistics //
  ht04ud_lr = new TH2F("ht04ud_lr", "dtud v.s. dtlr", 500, -8., 8., 500, -8., 8.);
  ht04ud_lr->GetXaxis()->SetTitle("tl-tr [ns]"); ht04ud_lr->GetYaxis()->SetTitle("tu-td [ns]");
  hx04ud_lr = new TH2F("hx04ud_lr", "Beam Spot @ F1 (given by dt)",
    300, -800., 800., 300, -800., 800.);
  hx04ud_lr->GetXaxis()->SetTitle("x [mm]"); hx04ud_lr->GetYaxis()->SetTitle("y [mm]");
  hlr_t02_t04 = new TH2F("hlr_t02_t04", "-dtlr_f1 v.s. dtlr_etf (daq alignment)",
    300, -7., 7., 300, -7., 7.);
  hlr_t02_t04->GetXaxis()->SetTitle("tl-tr_ETF [ns]");
  hlr_t02_t04->GetYaxis()->SetTitle("tr-tl_F1 [ns]");
  hx_f1_lr = new TH2F("hx_f1_lr", "x_f1 v.s. dtlr04", 52, -52., 52., 500, -1.7, 1.7);
  hx_f1_lr->GetXaxis()->SetTitle("x_F1 [mm]"); hx_f1_lr->GetYaxis()->SetTitle("tl-tr_F1 [ns]");
  hx_f1_ud = new TH2F("hx_f1_ud", "x_f1 v.s. dtud04", 52, -52., 52., 500, -4., 4.);
  hx_f1_ud->GetXaxis()->SetTitle("x_F1 [mm]"); hx_f1_ud->GetYaxis()->SetTitle("tu-td_F1 [ns]");
  hxf1_xETF = new TH2F("hxf1_xETF", "xETF v.s. xF1", 36, -45., 45., 52, -52., 52.);
  hxf1_xETF->GetXaxis()->SetTitle("x_ETF [mm]"); hxf1_xETF->GetYaxis()->SetTitle("x_F1 [mm]");
  hx_f1_tof1 = new TH2F("hx_f1_tof1", "tof1 v.s. x_f1", 52, -52., 52., 1000, 110., 180.);
  hx_f1_tof1->GetXaxis()->SetTitle("x_F1 [mm]"); hx_f1_tof1->GetYaxis()->SetTitle("tof1 [ns]");
  // fired distribution of the tofw and the t0arr //
  hx_f1 = new TH1F("hx_f1", "X Hit Position @ F1", 52, -52., 52.);
  hx_f1->GetXaxis()->SetTitle("x_F1 [mm]");
  htofwDis = new TH1F("htofwDis", "TOF Wall Fired Strip Distribution", 17, -1.5, 15.5);
  htofwDis->GetXaxis()->SetTitle("Strip Id");
  htofwDT = new TH2F("htofwDT", "TOF Wall TU-TD v.s. strip", 19, -1.5, 17.5, 500, -20., 20.);
  htofwDT->GetXaxis()->SetTitle("strip Id"); htofwDT->GetYaxis()->SetTitle("tu-td [ns]");
  hVETO = new TH2F("hVETO", "VETO HIT DISTRIBUTION", 5, -0.5, 4.5, 5, -0.5, 4.5);

  addh(hT0DtPvV, "daq"); addh(hT1DtPvV, "daq");
  addh(hT0LR, "daq"); addh(hT1LR, "daq"); addh(hT0LE); addh(hT0RE);
  addh(hT1LE); addh(hT1RE); addh(hT0L); addh(hT0R); addh(hT1L); addh(hT1R);
  addh(ht0arrMulti, "f1"); addh(htofwMulti, "tof"); addh(htof2, "tof");
  addh(htof1, "tof"); addh(htof1_v, "tof"); addh(htof1p_v, "tof");
  addh(htofwDis, "tof"); addh(htofwDT, "tof"); addh(hVETO, "beam");
  addh(hdE0, "ach"); addh(hdE1, "ach");
  addh(htof1_v_dE0, "pid"); addh(htof1_v_dE1, "pid");
  addh(htof1_dE0, "pid"); addh(htof1_dE1, "pid");
  addh(htof2_dE1, "pid"); addh(ht04ud_lr, "f1"); addh(hx04ud_lr, "f1");
  addh(hlr_t02_t04, "f1"); addh(hx_f1_lr, "f1"); addh(hx_f1_ud, "f1");
  addh(hxf1_xETF, "f1"); addh(hx_f1_tof1, "f1"); addh(hx_f1, "f1");


  // multiplicities and positions //
  const char xuv[3] = {'X', 'U', 'V'};
  char name[64]{}, title[128]{};
  // for dcs downstream of the dipole magnet
  for(int i = 0; i < 3; i++) for(int j = 0; j < 6; j++){
    const int dc = j / 2, x12 = j % 2;
    //--- multiplicity ---//
    sprintf(name, "hmdc%d%c%d", dc, xuv[i], x12);
    // only this manner of binning and title triggers automatic efficiency calculation
    // so do NOT bother change them to something else
    sprintf(title, "Length$(dc%d%c%d) -- Multiplicity", dc, xuv[i], x12);
    hmdc[i][j] = new TH1F(name, title, 13, -1.5, 11.5);
    addh(hmdc[i][j], "dcMul"); hmdc[i][j]->GetXaxis()->SetTitle("multiplicity");
    //--- fired position distribution ---//
    sprintf(name, "hpdc%d%c%d", dc, xuv[i], x12);
    sprintf(title, "Position Distribution -- dc%d%c%d", dc, xuv[i], x12);
    hpdc[i][j] = new TH1F(name, title, 83, -1.5, 81.5);
    addh(hpdc[i][j], "dcDis"); hpdc[i][j]->GetXaxis()->SetTitle("Wire Id");
    //--- drift time ---//
    sprintf(name, "htdc%d%c%d", dc, xuv[i], x12);
    sprintf(title, "Drift Time -- dc%d%c%d", dc, xuv[i], x12);
    htdc[i][j] = new TH1F(name, title, 500, -100., 400.);
    addh(htdc[i][j], "dcDT"); htdc[i][j]->GetXaxis()->SetTitle("drift time [ns]");
  } // end for over i and j
  hadcx = new TH1F("hadcx", "Angle in degree - postMag DCs", 500, 5., 30.);
  addh(hadcx, "trk"); hadcx->GetXaxis()->SetTitle("Angle X [degree]");
  // for dcs around the target
  const char xy[2] = {'X', 'Y'}, ud[2] = {'U', 'D'};
  for(int i = 0; i < 2; i++) for(int j = 0; j < 2; j++){
    for(int k = 0; k < 4; k++){
      //--- multiplicity ---//
      const int dc = k / 2, x12 = k % 2;
      sprintf(name, "hmdcTa%c%d%c%d", ud[i], dc, xy[j], x12);
      // only this manner of binning and title triggers automatic efficiency calculation
      // so do NOT bother change them to something else
      sprintf(title, "Length$(dcTa%c%d%c%d) -- Multiplicity", ud[i], dc, xy[j], x12);
      hmdcTa[i][j][k] = new TH1F(name, title, 13, -1.5, 11.5);
      addh(hmdcTa[i][j][k], "dcMul"); hmdcTa[i][j][k]->GetXaxis()->SetTitle("multiplicity");
      //--- fired position distribution ---//
      sprintf(name, "hpdcTa%c%d%c%d", ud[i], dc, xy[j], x12);
      sprintf(title, "Fired Wire Distribution -- dcTa%c%d%c%d", ud[i], dc, xy[j], x12);
      hpdcTa[i][j][k] = new TH1F(name, title, 18, -1.5, 16.5); // dcTaU
      addh(hpdcTa[i][j][k], "dcDis"); hpdcTa[i][j][k]->GetXaxis()->SetTitle("Wire Id");
      //--- drift time ---//
      sprintf(name, "htdcTa%c%d%c%d", ud[i], dc, xy[j], x12);
      sprintf(title, "Drift Time -- dcTa%c%d%c%d", ud[i], dc, xy[j], x12);
      htdcTa[i][j][k] = new TH1F(name, title, 500, -100., 400.);
      addh(htdcTa[i][j][k], "dcDT"); htdcTa[i][j][k]->GetXaxis()->SetTitle("drift time [ns]");
    } // end for over k
    // 2D versions //
    //--- fired position distribution ---//
    sprintf(name, "hp2dcTa%c%d", ud[i], j);
    sprintf(title, "Position Distribution Y vs X -- dcTa%c%d", ud[i], j);
    if(i) hp2dcTa[i][j] = new TH2F(name, title, 36, -90., 90., 36, -90., 90.); // dcTaD
    else hp2dcTa[i][j] = new TH2F(name, title, 36, -45., 45., 36, -45., 45.); // dcTaU
    addh(hp2dcTa[i][j], "pos");
    hp2dcTa[i][j]->GetXaxis()->SetTitle("X [mm]");
    hp2dcTa[i][j]->GetYaxis()->SetTitle("Y [mm]");
    //--- multiplicity ---//
    sprintf(name, "hm2dcTa%c%d", ud[i], j);
    sprintf(title, "Multiplicity Y vs X -- dcTa%c%d", ud[i], j);
    hm2dcTa[i][j] = new TH2F(name, title, 13, -1.5, 11.5, 13, -1.5, 11.5);
    addh(hm2dcTa[i][j], "dcMul");
    hm2dcTa[i][j]->GetXaxis()->SetTitle("Multiplicity - X1 || X2");
    hm2dcTa[i][j]->GetYaxis()->SetTitle("Multiplicity - Y1 || Y2");
    //--- angle ---//
    sprintf(name, "haTa%c%c", ud[i], xy[j]);
    sprintf(title, "Angle in Deg -- dcTa%c%c", ud[i], xy[j]);
    haTa[i][j] = new TH1F(name, title, 500, -5., 5.);
    sprintf(title, "Angle %c [degree]", xy[j]); haTa[i][j]->GetXaxis()->SetTitle(title);
    addh(haTa[i][j], "trk");
  } // end for over i and j
  // beam spot on the target //
  for(int i = 0; i < 2; i++){
    sprintf(name, "hTAp%c", ud[i]);
    sprintf(title, "Beam Spot On The Target by DCTa%c", ud[i]);
    hTAp[i] = new TH2F(name, title, 500, -90., 90., 500, -90., 90.); // dcTaD
    hTAp[i]->GetXaxis()->SetTitle("X [mm]");
    hTAp[i]->GetYaxis()->SetTitle("Y [mm]");
    addh(hTAp[i], "trk");
  } // end for over i
  //--- pid ---//
  for(int i = 0; i < 2; i++){
    // z
    sprintf(name, "hz%d", i);
    sprintf(title, "Z of MUSIC%d", i);
    hz[i] = new TH1F(name, title, 500, 0., 50.);
    hz[i]->GetXaxis()->SetTitle("Z");
    addh(hz[i], "pid");
  } // end for over i
  hpid0 = new TH2F("hpid0", "PID of Incident Beams", 500, 1.8, 2.5, 500, 0., 40.);
  addh(hpid0, "pid");
} // end member function CreateHistos

void ETFUserAnaTask::Configure(){
  fT1p = dynamic_cast<ETFTimer2 *>(fdm->GetDetector("tstop"));
  fdcTa[0][0] = dynamic_cast<ETFDCTa *>(fdm->GetDetector("dcTaU0"));
  fdcTa[0][1] = dynamic_cast<ETFDCTa *>(fdm->GetDetector("dcTaU1"));
  fdcTa[1][0] = dynamic_cast<ETFDCTa *>(fdm->GetDetector("dcTaD0"));
  fdcTa[1][1] = dynamic_cast<ETFDCTa *>(fdm->GetDetector("dcTaD1"));
  fdc[0] = dynamic_cast<ETFDC *>(fdm->GetDetector("dc0"));
  fdc[1] = dynamic_cast<ETFDC *>(fdm->GetDetector("dc1"));
  fdc[2] = dynamic_cast<ETFDC *>(fdm->GetDetector("dc2"));
  fT1v = dynamic_cast<ETFTimer2 *>(fdm->GetDetector("tstop0_v"));
  fVETO = dynamic_cast<ETFTimer4 *>(fdm->GetDetector("veto_v"));
  fT0p = dynamic_cast<ETFTimer2 *>(fdm->GetDetector("tstart2"));
  fT0v = dynamic_cast<ETFTimer2 *>(fdm->GetDetector("tstart2_v"));
  fTOFW = dynamic_cast<ETFTOFWall *>(fdm->GetDetector("tofwB"));
  if(!fTOFW) fTOFW = dynamic_cast<ETFTOFWall *>(fdm->GetDetector("tofw"));
  fT0F1 = dynamic_cast<ETFTimer4 *>(fdm->GetDetector("tstart4"));
  fT0Arr = dynamic_cast<ETFT0Array *>(fdm->GetDetector("t0arr"));
  fMU[0] = dynamic_cast<ETFMUSIC *>(fdm->GetDetector("mu0"));
  fMU[1] = dynamic_cast<ETFMUSIC *>(fdm->GetDetector("mu1"));
} // end member function Configure

void ETFUserAnaTask::Exec(){
  Reset(); // reset related member data variables

  //----- T0 Array @F1 related statistics -----//
  double dtlr04 = -999., dtud04 = -999.;
  static const double xeff = 94.87; // =VEFF/2, in mm/ns
  if(fT0F1){
    // tstart4 @F1
    ETFPlaHit *hlr04 = fT0F1->GetHit(0);
    ETFPlaHit *hud = fT0F1->GetHit(1);
    // stripId of lr should be 0
    if((hlr04 && hlr04->stripId() == 1) || (hud && hud->stripId() == 0)) swap(hlr04, hud);
    if(hlr04 && hlr04->fired() == 4) dtlr04 = hlr04->dt();
    if(hud && hud->fired() == 4) dtud04 = hud->dt();
    if(-999. != dtlr04 && -999. != dtud04){
      ht04ud_lr->Fill(dtlr04, dtud04);
      hx04ud_lr->Fill(dtlr04*xeff, dtud04*xeff);
    } // end if
    // tstart2 to ETF
    if(fT0p){
      double dtlr02 = -999.;
      ETFPlaHit *hlr = fT0p->GetHit(0);
      if(hlr && hlr->fired() == 4) dtlr02 = hlr->dt();
      if(-999. != dtlr04 && -999. != dtlr02) hlr_t02_t04->Fill(-dtlr04, dtlr02);
    } // end if(fT0p)
    if(fT0F1->GetFiredStatus() == 4) t0_f1 = fT0F1->GetTime();
  } // end if
  if(fT0Arr){
    // if(t0_f1 > 0.) fT0Arr->SetT0(t0_f1);
    ETFPlaHit *tm = fT0Arr->OptimalHit();
    if(tm){
      x_f1 = tm->x();
      t0arr_m = tm->time();
    } // end if
    if(-999. != x_f1){
      if(-999. != dtlr04) hx_f1_lr->Fill(x_f1, dtlr04);
      if(-999. != dtud04) hx_f1_ud->Fill(x_f1, dtud04);
    } // end if
    ht0arrMulti->Fill(fT0Arr->multi());
  } // end if(fT0Arr)

  // PXI v.s. VME alignment check, plus some t0,1 statistics //
  double tp0[2], tp1[2]; // T0,1[L,R]
  for(int i = 2; i--;) tp0[i] = tp1[i] = -999.;
  // assign the data //
  if(fPXISections.size() && fPXISections[0]->UnpackGood()){
    if(fT0p){ tp0[0] = fT0p->timeL(); tp0[1] = fT0p->timeR(); }
    if(fT1p){ tp1[0] = fT1p->timeL(); tp1[1] = fT1p->timeR(); }
    if(fT0p && fT1p && fT0p->fired() == 4 && fT1p->fired() == 4)
      tof1 = fT1p->time() - fT0p->time();
    // never use tof by subtraction of two daqs - 25ns resolution, very bad //
    // if(fT0F1 && fT0F1->fired() == 4 && fT1p && fT1p->fired() == 4)
    //   tof1 = fT1p->time() - fT0F1->time() + 1771.5;
    htof1->Fill(tof1);
  } // end if
  double t0_v = -999., t1_v = -999.; // tstart and tstop measured by VME
  if(fT0v){ // t0 daqed by vme: time v.s. edges
    TClonesArray *c = fT0v->GetHits();
    if(c->GetEntriesFast()){
      ETFPlaHit *h = dynamic_cast<ETFPlaHit *>(c->At(0));
      const ETFHPTDCCh *tl = h->GetTChannel(0); // up, or left
      const ETFHPTDCCh *tr = h->GetTChannel(2); // down, or right
      if(tl && tr){
        const int nlL = tl->nl(), nlR = tr->nl(), nl = nlL <= nlR ? nlL : nlR;
        double tm[2] = {-999., -999.}; // the edges that match PXI daq
        for(int i = 0; i < nlL; i++){
          const double t = tl->GetLeadingEdge(i);
          hT0LE->Fill(i, t);
          hT0L->Fill(t);
          // vme mtdc takes a very long match window (46us) and trigger latency (24us)
          // [-23000, trig, 23000], so t-22600 < 500 selects trigger-generating hit
          // taken as the trigger-generating hit
          if(-999. == tm[0]){
            if(ETFMath::Within(t, 22200., 25000.)) tm[0] = t;
          } // end if
        } // end for over i
        for(int i = 0; i < nlR; i++){
          const double t = tr->GetLeadingEdge(i);
          hT0RE->Fill(i, t);
          hT0R->Fill(t);
          // taken as the trigger-generating hit
          if(-999. == tm[1]){
            if(ETFMath::Within(t, 22200., 25000.)) tm[1] = t;
          } // end if
        } // end for over i
        if(-999. != tm[0] && -999. != tm[1]){
          t0_v = (tm[0] + tm[1]) * 0.5 - fT0v->GetDelay();
          if(-999. != tp0[0] && -999. != tp0[1])
            hT0DtPvV->Fill(tp0[0] - tp0[1], tm[0] - tm[1]);
        } // end if
        for(int i = 0 ; i < nl; i++)
          hT0LR->Fill(tl->GetLeadingEdge(i), tr->GetLeadingEdge(i));
      } // end if
    } // end if
  } // end T0
  if(fT1v){ // time v.s. edges - tstop
    TClonesArray *c = fT1v->GetHits();
    if(c->GetEntriesFast()){
      ETFPlaHit *h = dynamic_cast<ETFPlaHit *>(c->At(0));
      const ETFHPTDCCh *tl = h->GetTChannel(0); // up, or left
      const ETFHPTDCCh *tr = h->GetTChannel(2); // down, or right
      if(tl && tr){
        const int nlL = tl->nl(), nlR = tr->nl(), nl = nlL <= nlR ? nlL : nlR;
        double tm[2] = {-999., -999.}; // the vme edges that match PXI daq
        for(int i = 0; i < nlL; i++){
          const double t = tl->GetLeadingEdge(i);
          hT1LE->Fill(i, t);
          hT1L->Fill(t);
          // taken as the trigger-generating hit
          if(-999. == tm[0]){
            if(ETFMath::Within(t, 22200., 25000.)) tm[0] = t;
          } // end if
        } // end for over i
        for(int i = 0; i < nlR; i++){
          const double t = tr->GetLeadingEdge(i);
          hT1RE->Fill(i, t);
          hT1R->Fill(t);
          // taken as the trigger-generating hit
          if(-999. == tm[1]){
            if(ETFMath::Within(t, 22200., 25000.)) tm[1] = t;
          } // end if
        } // end for over i
        if(-999. != tm[0] && -999. != tm[1]){
          t1_v = (tm[0] + tm[1]) * 0.5 - fT1v->GetDelay();
          if(-999. != tp1[0] && -999. != tp1[1])
            hT1DtPvV->Fill(tp1[0] - tp1[1], tm[0] - tm[1]);
        } // end if
        for(int i = 0 ; i < nl; i++)
          hT1LR->Fill(tl->GetLeadingEdge(i), tr->GetLeadingEdge(i));
      } // end if
    } // end if
  } // end T1
  if(-999. != t0_v && -999. != t1_v){
    tof1_v = t1_v - t0_v;
    htof1_v->Fill(tof1_v);
    if(-999. != tof1) htof1p_v->Fill(tof1, tof1_v);
  } // end if
  // VETO statistics
  if(fVETO){
    const ETFHPTDCCh *t[4]{};
    TClonesArray *cv = fVETO->GetHits();
    const int np = cv->GetEntriesFast();
    for(int j = 0; j < np; j++){
      if(j > 1) continue; // not possible
      ETFPlaHit *h = dynamic_cast<ETFPlaHit *>(cv->At(j));
      t[j*2] = h->tc(0); // UV
      t[j*2+1] = h->tc(2); // DV
    } // end for over j
    // beam position //
    // off-target
    bool fired[4]{};
    for(int j = 4; j--;) if(t[j] && t[j]->nl()) fired[j] = true;
    if(fired[0]) hVETO->Fill(2, 3); // up
    else if(fired[1]) hVETO->Fill(1, 2); // left (x<0)
    else if(fired[2]) hVETO->Fill(2, 1); // down
    else if(fired[3]) hVETO->Fill(3, 2); // right (x>0)
    // in the bull's eye
    else hVETO->Fill(2, 2);
  } // end if(veto)

  if(fTOFW){
    const int mul = fTOFW->multi();
    htofwMulti->Fill(mul);
    for(TObject *b : *fTOFW->hits()){
      ETFPlaHit *h = dynamic_cast<ETFPlaHit *>(b);
      const int id = h->stripId();
      if(h && h->fired() == 4){
        htofwDis->Fill(id);
        htofwDT->Fill(id, h->dt());
        if(1 == mul) if(fT1p && fT1p->fired() == 4){
          const double t1 = fT1p->time();
          tof2 = h->GetTime() - t1 + fPID->tof2_ofst();
          beta = 12000. / (tof2*ETFMath::c0());
          if(beta > 0. && beta < 1.) bg = ETFMath::BetaGamma(beta);
          htof2->Fill(id, tof2);
        } // end if t1
      } // end if
    } // end for over tofw strips
  } // end if TOFW

  //----- tracking related -----//
  if(fAnalyzer->AnaLevel() >= 1){ // tracking is on, only used in offline mode
    // 1. get the trackers //
    ETFTaTracker *preTaTracker = ftm->GetPreTaTracker();
    ETFTaTracker *postTaTracker = ftm->GetPostTaTracker();
    ETFPostMagTracker *postMagTracker = ftm->GetPostMagTracker();
    // 2. get the tracks //
    ETFVTrack *preTaTrk[2]{0}, *postTaTrk[2]{0}; // [x-y]
    ETFVTrack *postMagTrk = nullptr;
    ETF3DTrack *postMag3DTrk = nullptr;
    if(preTaTracker){
      if(preTaTracker->GetNofTracks(TrackType::kX) == 1)
        preTaTrk[0] = preTaTracker->GetTrack(0, TrackType::kX);
      if(preTaTracker->GetNofTracks(TrackType::kY) == 1)
        preTaTrk[1] = preTaTracker->GetTrack(0, TrackType::kY);
    } // end if
    if(postTaTracker){
      if(postTaTracker->GetNofTracks(TrackType::kX) == 1)
        postTaTrk[0] = postTaTracker->GetTrack(0, TrackType::kX);
      if(postTaTracker->GetNofTracks(TrackType::kY) == 1)
        postTaTrk[1] = postTaTracker->GetTrack(0, TrackType::kY);
    } // end if
    if(postMagTracker){
      if(postMagTracker->GetNofTracks(TrackType::kX) == 1)
        postMagTrk = postMagTracker->GetTrack(0, TrackType::kX);
      if(postMagTracker->GetNof3DTracks() == 1)
        postMag3DTrk = postMagTracker->Get3DTrack(0);
    } // end if
    // 3. extract information
    for(int i = 0; i < 2; i++){ // loop over x and y
      if(preTaTrk[i]){
        k0[i] = preTaTrk[i]->GetK();
        b0[i] = preTaTrk[i]->GetB();
        pT1[i] = preTaTrk[i]->GetHitPos("tstop");
        pTAU[i] = preTaTrk[i]->GetHitPos("TA");
        pdcTaU[0][i] = preTaTrk[i]->GetHitPos("dcTaU0");
        pdcTaU[1][i] = preTaTrk[i]->GetHitPos("dcTaU1");
        hTAp[0]->Fill(pTAU[0], pTAU[1]);
        haTa[0][i]->Fill(preTaTrk[i]->ang());
      } // end if
      if(postTaTrk[i]){
        k1[i] = postTaTrk[i]->GetK();
        b1[i] = postTaTrk[i]->GetB();
        pTAD[i] = postTaTrk[i]->GetHitPos("TA");
        pdcTaD[0][i] = postTaTrk[i]->GetHitPos("dcTaD0");
        pdcTaD[1][i] = postTaTrk[i]->GetHitPos("dcTaD1");
        hTAp[1]->Fill(pTAD[0], pTAD[1]);
        haTa[1][i]->Fill(postTaTrk[i]->ang());
      } // end if
    } // end for over x and y
    if(postMagTrk){
      k2[0] = postMagTrk->GetK();
      b2[0] = postMagTrk->GetB();
      pdc[0][0] = postMagTrk->GetHitPos("dc0");
      pdc[1][0] = postMagTrk->GetHitPos("dc1");
      pdc[2][0] = postMagTrk->GetHitPos("dc2");
      ptofw[0] = postMagTrk->GetHitPos("tofwB");
      hadcx->Fill(postMagTrk->ang());
    } // end if
    if(postMag3DTrk){
      k2[1] = postMag3DTrk->GetK2();
      b2[1] = postMag3DTrk->GetB2();
      pdc[0][1] = postMag3DTrk->GetHitPosY("dc0");
      pdc[1][1] = postMag3DTrk->GetHitPosY("dc1");
      pdc[2][1] = postMag3DTrk->GetHitPosY("dc2");
      ptofw[1] = postMag3DTrk->GetHitPosY("tofwB");
    } // end if
  } // end if tracking is on
  
  // DCTa statistics //
  int mdcTa[2][2][4]{}; // [U-D][X-Y][dcTaU(D)0-x0-1--dcTaU(D)1] - multi
  TClonesArray *cl;
  for(int l = 0; l < 2; l++){ // loop over U-D
    for(int i = 0; i < 2; i++){ // loop over dcTaU(D)0-1
      for(int j = 0; j < 2; j++){ // loop over x-y
        pdcTaWire[l][i][j] = fdcTa[l][i]->WireHitPos(j);
        for(int k = 0; k < 2; k++){ // loop over x0-1
          cl = fdcTa[l][i]->GetHitsArray(DCType(j), k);
          const int nh = mdcTa[l][j][i*2+k] = cl->GetEntriesFast();
          hmdcTa[l][j][i*2+k]->Fill(nh);
          for(int m = 0; m < nh; m++){
            ETFDCHit *h = dynamic_cast<ETFDCHit *>(cl->At(m));
            if(1 == nh) hpdcTa[l][j][i*2+k]->Fill(h->GetWireId());
            if(h->fired() >= 2) htdcTa[l][j][i*2+k]->Fill(h->GetT());
          } // end for over m
        } // end for over k
      } // end for over j
      hp2dcTa[l][i]->Fill(pdcTaWire[l][i][0], pdcTaWire[l][i][1]);
      for(int j = 2; j--;) hm2dcTa[l][i]->Fill(mdcTa[l][0][i*2+j], mdcTa[l][1][i*2+j]);
    } // end for over i
  } // end for over l
  // postMag dc statistics //
  for(int i = 0; i < 3; i++){ // X-U-V
    for(int j = 0; j < 6; j++){ // dc0-x0-1--dc1--dc2
      const int dc = j / 2, x12 = j % 2;
      cl = fdc[dc]->GetHitsArray(DCType(i), x12);
      const int nh = cl->GetEntriesFast();
      hmdc[i][j]->Fill(nh);
      for(int k = 0; k < nh; k++){
        ETFDCHit *h = dynamic_cast<ETFDCHit *>(cl->At(k));
        if(1 == nh) hpdc[i][j]->Fill(h->GetWireId());
        if(h->fired() >= 2) htdc[i][j]->Fill(h->GetT());
      } // end for over k
    } // end for over j
  } // end for over i

  static const double LL = ETFVPID::L / ETFMath::c0();
  if(-999. != tof1){
    beta0 = LL / tof1;
    brho0 = fBrho0;
    if(beta0 > 0. && beta0 < 1.) bg0 = ETFMath::BetaGamma(beta0);
    if(bg0 > 0. && x_f1 != -999. && pdcTaWire[0][1][0] != -999.){
      // 30Dec, Ni-58 setting, Brho = 4.9896745, brho0*CC=1.6058799
      // CC = 0.32184060506023993,  = e0/(u0*c0*1E6), in SI unit
      aoz0 = brho0*ETFMath::CC()/bg0;
      // the guessed one - more separated, yet aoz distored --- :) current preferred one
      aoz0 = aoz0*(1.+0.000194047*pdcTaWire[0][1][0]-0.0011855*x_f1); // correction
      aoz0 = aoz0+0.19894667; // calibration
      // the fitted one - more physical, yet less separated -- :( less favored
      // aoz0 = aoz0*(1.+0.000207*pdcTaWire[0][1][0]-0.00816106*x_f1); // correctionex
      // aoz0 = aoz0*1.03974-0.130472; // calibration
    } // end if
  } // end if
  //----- correct tof1 with what we've all got -----//
  if(-999. != x_f1){
    if(-999. != pdcTaWire[0][1][0]) hxf1_xETF->Fill(pdcTaWire[0][1][0], x_f1);
    if(-999. != tof1) hx_f1_tof1->Fill(x_f1, tof1);
    hx_f1->Fill(x_f1);
  } // end if

  //----- Delta E and Z related -----//
  for(int i = 0; i < 2; i++) if(fMU[i]) dE[i] = fMU[i]->dE();
  if(dE[0] > 0.){
    hdE0->Fill(dE[0]);
    if(pdcTaWire[0][1][1] != -999. && tof1 > 0.){
      // B-B formula with Y correction //
      // const double de = dE[0]*exp(0.0017*(pdcTaWire[0][1][1]+30.)); // pos correction
      // 2023dec
      // z0 = sqrt(de/(0.00228*pow(tof1-121.43353,2)+6.973)); // beta correction
      // z0 = z0*63.3748-1.22128; // calibration
      // 2024jun
      z0 = sqrt(dE[0]*exp(0.0025*(pdcTaWire[0][1][1]+30.))/(0.0023*pow(tof1-121.,2)+7.))*53.9046-0.480279;
      hz[0]->Fill(z0);
      // if(aoz0 > 0.){
      //   z0c = z0Straight(aoz0, z0); // make z0 straight
      //   aoz0c = aoz0Straight(aoz0, z0); // give aoz0 physical meaning
      // } // end if
    } // end if
    htof1_v_dE0->Fill(tof1_v, dE[0]);
    htof1_v_dE1->Fill(tof1_v, dE[1]);
    htof1_dE0->Fill(tof1, dE[0]);
    htof1_dE1->Fill(tof1, dE[1]);
  } // end if
  if(dE[1] > 0.){
    hdE1->Fill(dE[1]);
    htof2_dE1->Fill(tof2, dE[1]);
    // z calibration //
    if(beta > 0. && beta < 1. && pdcTaWire[1][1][1] != -999.){
      const double de = dE[1]*exp(0.00035*(pdcTaWire[1][1][1]+40.)); // pos correction
      z1 = sqrt(de)+0.63*beta-0.2*beta*de; // transform to a z-linear variable
      z1 = z1*33.3786-18.2236; // calibrate to Z
      // FIXME FIXME - provisional //
      z1 = z1*pow(beta0,0.3)*1.30885-1.08661; // only valid for NOTA Exp - temporary
      hz[1]->Fill(z1);
    } // end if
  } // end if
  if(fAnalyzer->AnaLevel() < 2){
    if(tof1 > 0.) fPID->SetTOF1(tof1);
    if(tof2 > 0.) fPID->SetTOF2(tof2);
  } // end if
  if(-999. != aoz0 && -999. != z0) hpid0->Fill(aoz0, z0);

  //----- pileup related -----//
  if(fT1v){
    if(fT1v->fired() > 1){
      ETFPlaHit *h = dynamic_cast<ETFPlaHit *>(fT1v->GetHits()->At(0));
      const double trigLB = 22750. - 1000., trigHB = 22750. + 1000.;
      // count the particles that echo in both L and R (fired=4) in fPileUpN
      // and determine min dt [fabs(tl-tr)] of the particles above which hit
      // the det within (trigLB, trigHB) in fPileUpT
      h->PileUpAna(fPileUpN, fPileUpT, trigLB, trigHB);
    } // end if(fT1v->fired() > 1)
  } // end if(fT1v)

  // post-ta pid
  if(fAnalyzer->AnaLevel() >= 2 && fPID->tof2() > 0.){
    tof2 = fPID->tof2(); // refined tof2
    beta = fPID->beta(); // and beta
    rho = fPID->rho();
    aoz = fPID->aoz();
  } // end if

  et->Fill();
} // end member function Exec

void ETFUserAnaTask::WriteHistos(){
  ETFTask::WriteHistos();
  fTree->AddFriend(et);
  et->Write("", TObject::kOverwrite);
} // end member function WriteHistos

void ETFUserAnaTask::Reset(){
  x_f1 = t0_f1 = t0arr_m = -999.;
  tof1 = tof1_v = tof2 = -999.;
  z0 = z1 = -999.;
  fPileUpN = -1; fPileUpT = -999.;
  for(int i = 0; i < 2; i++){
    k0[i] = b0[i] = -999.; // preTa
    k1[i] = b1[i] = -999.; // postTa
    k2[i] = b2[i] = -999.; // postTa
  } // end for over i
  aoz0 = aoz = -999.;
  brho0 = bg0 = beta0 = beta = bg = -1.;
  rho = -999.;
  aoz0c = -999.; z0c = -1.;

  for(int i = 0; i < 2; i++){
    dE[i] = -999.;
    pT1[i] = pTAU[i] = pTAD[i] = -999.; // U(D): from tracking of dcTaU(D) array
    for(int j = 0; j < 2; j++){
      pdcTaU[i][j] = pdcTaD[i][j] = -999.;
      pdcTaWire[0][i][j] = pdcTaWire[1][i][j] = -999.;
    } // end for over j
    pdc[0][i] = pdc[1][i] = pdc[2][i] = ptofw[i] = -999.;
  } // end for over i
} // end member function Reset


// the piecewise fitted function - returns the fitted value
// while the intended value is returned in Y
// ii specifies the piece of x range
double piecefit(int ii, double x, double y, double &Y);
// x,y: aoz0 and z0
double ETFUserAnaTask::z0Straight(double x, double y){
  double y0, y1, Y0, Y1; // y0,1: the fitted values; Y0,1: the intended values
  const int yi = int(y+0.5); // the main side
  y0 = piecefit(yi, x, y, Y0);
  // if(-1 == y0) return y; // indiscernible - do nothing
  // if(-2 == y0) return y; // x unavaliable - do nothing
  const int yi1 = y > y0 ? yi + 1 : yi - 1; // the other side
  y1 = piecefit(yi1, x, y, Y1);
  // either side unavailable - suppose the other side parallel with this side
  // if(y1 < 0.){ y1 = y0 + yi1 - yi; Y1 = yi1; }
  // if(y0 < 0.){ y0 = y1 - yi1 + yi; Y0 = yi;  }

  // now we do the linear interpolation
  return lii(y0, Y0, y1, Y1, y);
} // end member function z0Straight

// \retval -1: z out of border; -2: x out of border
double piecefit(int ii, double x, double y, double &Y){
  // depict the pid spectrum //
  static const int NZ = 19, NI = 4; // NI: # of isotopes per Z
  static const int ZMAX = 37, ZMIN = ZMAX - NZ + 1;
  static const double p[NZ][NI][2][2] = { // [z][isotope][2ends][(x,y)]
    { // z = 37
      {{-1., 37.}, {-1., 37.}}, // Rb74
      {{-1., 37.}, {-1., 37.}}, // Rb75
      {{2.04287, 37.4409}, {2.0776, 36.9015}}, // Rb76
      {{2.07508, 37.3061}, {2.11686, 36.9}} // Rb77
    },
    { // z = 36
      // {{-1., 36.}, {-1., 36.}}, // Kr72
      {{2.01619, 36.32283}, {2.04589, 36.0923}}, // Kr73
      {{2.04438, 36.22}, {2.08012, 35.84}}, // Kr74
      {{2.07508, 36.2}, {2.11686, 35.75}}, // Kr75
      {{2.14328, 36.337}, {2.18655, 36.}} // Kr78
    },
    { // z = 35
      {{-1., 35.}, {-1., 35.}}, // Br70
      {{2.01367, 35.1819}, {2.04488, 34.8448}}, // Br71
      {{2.04438, 35.13}, {2.08213, 34.87}}, // Br72
      {{2.08012, 35.1}, {2.11837, 34.8}} // Br73
    },
    { // z = 34
      {{-1., 34.}, {-1., 34.}}, // Se68
      {{2.0182, 34.1018}, {2.04941, 33.852}}, // Se69
      {{2.04841, 34.1642}, {2.08616, 33.7896}}, // Se70
      {{2.08263, 34.22}, {2.11787, 33.78}} // Se71
    },
    { // z = 33
      {{-1., 33.}, {-1., 33.}}, // As66
      {{2.0172, 33.134}, {2.05143, 32.9155}}, // As67
      {{2.04891, 33.23}, {2.08867, 32.77}}, // As68
      {{2.08465, 33.10}, {2.11938, 32.90}} // As69
    },
    { // z = 32
      {{-1., 32.}, {-1., 32.}}, // Ge64
      {{2.01518, 32.19}, {2.05193, 31.81}}, // Ge65
      {{2.04841, 32.1662}, {2.08918, 31.7291}}, // Ge66
      {{2.08666, 32.1}, {2.12492, 31.7}} // Ge67
    },
    { // z = 31
      {{-1., 31.}, {-1., 31.}}, // Ga62
      {{2.01317, 31.26}, {2.05243, 30.74}}, // Ga63
      {{2.04891, 31.15}, {2.09018, 30.84}}, // Ga64
      {{2.08817, 31.1048}, {2.1224, 30.7614}} // Ga65
    },
    { // z = 30
      {{-1., 30.}, {-1., 30.}}, // Zn60
      {{2.01267, 30.2931}, {2.05042, 29.8248}}, // Zn61
      {{2.04941, 30.15}, {2.0922, 29.75}}, // Zn62
      {{2.08817, 30.1}, {2.1219, 29.9}} // Zn63
    },
    { // z = 29
      {{-1., 29.}, {-1., 29.}}, // Cu58
      {{2.01468, 29.138}, {2.05243, 28.7634}}, // Cu59
      {{2.05193, 29.138}, {2.09522, 28.7634}}, // Cu60
      {{2.09169, 29.1068}, {2.12643, 28.7634}} // Cu61
    },

    { // z = 28
      {{-1., 28.}, {-1., 28.}}, // Ni56
      {{2.01279, 28.2058}, {2.05557, 27.7975}}, // Ni57
      {{2.05416, 28.1183}, {2.09271, 27.7684}}, // Ni58
      {{2.09271, 28.1}, {2.12515, 27.9}} // Ni59
    },
    { // z = 27
      {{1.97518, 27.0977}, {2.01291, 26.8108}}, // Co54
      {{2.01291, 27.15}, {2.05538, 26.85}}, // Co55
      {{2.05395, 27.15}, {2.09356, 26.85}}, // Co56
      {{2.09165, 27.13}, {2.12792, 26.85}} // Co57
    },
    { // z = 26
      {{1.97808, 26.1301}, {2.0153, 25.8637}}, // Fe52
      {{2.01291, 26.2}, {2.05682, 25.8}}, // Fe53
      {{2.05347, 26.05}, {2.0969, 25.85}}, // Fe54
      {{2.0969, 26.1}, {2.12935, 25.9}} // Fe55
    },
    { // z = 25
      {{1.97903, 25.2126}, {2.01053, 24.8279}}, // Mn50
      {{2.01291, 25.2}, {2.05586, 24.8}}, // Mn51
      {{2.05538, 25.15}, {2.09928, 24.75}}, // Mn52
      {{2.10024, 25.1238}, {2.1365, 24.8126}} // Mn53
    },
    { // z = 24
      {{1.97501, 24.236}, {2.01207, 23.8216}}, // Cr48
      {{2.01772, 24.1176}, {2.05603, 23.8512}}, // Cr49
      {{2.06169, 24.088}, {2.09937, 23.94}}, // Cr50
      {{2.1, 24.1472}, {2.13643, 23.9104}} // Cr51
    },
    { // z = 23
      {{1.97187, 23.1409}, {2.01332, 22.845}}, // V46
      {{2.01646, 23.1409}, {2.05603, 22.8154}}, // V47
      {{2.06106, 23.0521}, {2.1, 22.8154}}, // V48
      {{2.10754, 23.1113}, {2.13769, 22.845}} // V49
    },
    { // z = 22
      {{1.97689, 22.0755}, {2.01521, 21.8979}}, // Ti44
      {{2.0196, 22.1347}, {2.05666, 21.8979}}, // Ti45
      {{2.06294, 22.1643}, {2.10251, 21.8979}}, // Ti46
      {{2.10942, 22.0755}, {2.14459, 21.8979}} // Ti47
    },
    { // z = 21
      {{1.97627, 20.9805}, {2.01081, 20.9509}}, // Sc42
      {{2.0196, 21.0988}, {2.05917, 20.8325}}, // Sc43
      {{2.06545, 21.1284}, {2.10691, 20.9213}}, // Sc44
      {{2.11382, 21.0396}, {2.14962, 20.9805}} // Sc45
    },
    { // z = 20
      {{1.97564, 19.9}, {2.01144, 19.8854}}, // Ca40
      {{2.02149, 20.0334}, {2.06043, 19.8262}}, // Ca41
      {{2.07048, 20.1518}, {2.10628, 19.8854}}, // Ca42
      {{2.11947, 20.0926}, {2.1603, 19.8558}} // Ca43
    },
    { // z = 19
      {{1.97742, 19.0231}, {2.01222, 19.0855}}, // K38
      {{2.01898, 19.0231}, {2.06208, 18.9607}}, // K39
      {{2.07351, 19.0543}, {2.11662, 18.9607}}, // K40
      {{2.12493, 18.9607}, {2.15557, 18.9919}} // K41
    }
  };

  if(!within(ii, ZMIN, ZMAX)){ Y = ii; return ii; } // out of Z zone - do nothing
  const int j = ZMAX - ii;
  // model the fitting function //
  // 1. pinpoint the isotope
  double ym; // the optimal one
  // if the 2 isotopes overlap a little, we gotta make a choice
  int ia[2] = {-1, -1}; // isotope id - primary, auxiliary
  for(int i = 0; i < NI; i++){
    if(-1 != p[j][i][0][0] && -1 != p[j][i][1][0] &&
        within(x, p[j][i][0][0], p[j][i][1][0])){
      if(-1 == ia[0]) ia[0] = i;
      else if(-1 == ia[1]) ia[1] = i;
      else break;
    } // end if
  } // end for over i
  if(-1 == ia[1]){
    if(-1 == ia[0]){ Y = ii; return ii; } // out of x range - do nothing
    ym = lii(p[j][ia[0]][0], p[j][ia[0]][1], x);
  } // end if
  else{ // select the isotope whose yy is closer to y
    double yy[2] = {-999., -999.};
    for(int i = 0; i < 2; i++) if(-1 != ia[i])
      yy[i] = lii(p[j][ia[i]][0], p[j][ia[i]][1], x);
    if(fabs(y-yy[0]) < fabs(y-yy[1])) ym = yy[0];
    else ym = yy[1];
  } // end else
  // 2. calculate the intended Y by correcting out the linear increment
  Y = ii;

  return ym;
} // end function piecefit

// give aoz0 physical meaning
// x,y: aoz0, z0
double ETFUserAnaTask::aoz0Straight(double x, double y){
  // depict the pid spectrum //
  static const int NZ = 22, ZMAX = 37;
  static const double p[NZ][2] = { // [z][p0-p1]
		{0.386077, 0.82249}, // z=37
		{0.442612, 0.794722}, // z=36
		{0.330241, 0.849819}, // z=35
		{0.345807, 0.843227}, // z=34
		{0.321307, 0.856027}, // z=33
		{0.265542, 0.883083}, // z=32
		{0.22359, 0.904365}, // z=31
		{0.212571, 0.91068}, // z=30
		{0.20509, 0.915439}, // z=29
		{0.153493, 0.941874}, // z=28
		{0.123823, 0.957624}, // z=27
		{0.100303, 0.970141}, // z=26
		{0.0733221, 0.984701}, // z=25
		{0.0360971, 1.00383}, // z=24
		{0.0173611, 1.01406}, // z=23
		{-0.0455366, 1.04664}, // z=22
		{0.00785906, 1.02162}, // z=21
		{-0.00820958, 1.03089}, // z=20
		{-0.0309923, 1.04419}, // z=19
		{0.0114466, 1.02479}, // z=18
		{-0.0269599, 1.04442}, // z=17
		{0.0688575, 0.997961} // z=16
  };
  const int ii = ZMAX - int(y+0.5);
  if(ii < 0 || ii >= NZ) return x;
  const double res = p[ii][0] + p[ii][1]*x;
  return res;
} // end member function aoz0Straight
